######################################################## 
##  Program Name  :/projects/bsi/bioinf_int/s106381.borawork/beagle_test/rpgm/compare_bglprob_geno.R
##  Study Title   : BORA
##  Programmer    : MEM
##  Lead :          de Andrade/Sicotte
##  Creation eDate: Monday, 21 February 2011 04:44 PM CST
##  Function      : Compare BEAGLE imputed genotype probabilities to actual genotype
######################################################## 
indir <- "/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr10"

##Get gprobs
setwd(indir)

##Get whole chrom run
whole <- read.table("whole_default/beagle_out_chr10.geno.bgl.gz.gprobs", header=T, as.is=T)
whole[,seq(4,ncol(whole),by=3)] <- 2*whole[,seq(4,ncol(whole),by=3)]+whole[,seq(5,ncol(whole),by=3)]
whole <- whole[,c(1:3, seq(4,ncol(whole),by=3))]

##Get whole chrom run with different seed
whole2 <- read.table(gzfile("whole_random2/beagle_out_chr10.geno.bgl.gz.gprobs"), header=T, as.is=T)
whole2[,seq(4,ncol(whole2),by=3)] <- 2*whole2[,seq(4,ncol(whole2),by=3)]+whole2[,seq(5,ncol(whole2),by=3)]
whole2 <- whole2[,c(1:3, seq(4,ncol(whole2),by=3))]

##Read in snpinfo
map <- read.table("/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr10/chr10_map.txt", header=T, as.is=T)
names(map) <- c("rsid", "bp")
map$chr <- 10
r2 <- read.table("whole_default/beagle_out_chr10.geno.bgl.gz.r2", header=F, as.is=T)
map <- map[map$rsid %in% whole$marker,]
map$r2 <- r2$V2[match(map$rsid, r2$V1)]
map[1:5,]
map$alleleA <- whole$alleleA[match(map$rsid, whole$marker)]
map$alleleB <- whole$alleleB[match(map$rsid, whole$marker)]
map <- map[match(whole$marker, map$rsid),]

#gprob <- system("find . -name \\*whole_gprobs",intern=T)
#gprob <- gprob[-grep("(seed|new)", gprob)]
#gprob
#gprob <- gprob[c(1,2,4)]
gprob <-c("./1mb/whole_gprobs","./2mb/whole_gprobs","./4mb/whole_gprobs","./halfmb/whole_gprobs","./halfquatermb/whole_gprobs","./quatermb/whole_gprobs","./quaterquatermb/whole_gprobs")
for (g in gprob) {
  print(g)
  sp <- gsub(".*[[:punct:]]", "", gsub("mb.*", "mb", g))
  temp <- read.table(g, header=T, as.is=T)
  temp[,seq(4,ncol(temp),by=3)] <- 2*temp[,seq(4,ncol(temp),by=3)]+temp[,seq(5,ncol(temp),by=3)]
  temp <- temp[,c(1:3, seq(4,ncol(temp),by=3))]
  temp <- temp[,names(whole)]
  temp <- temp[match(whole$marker, temp$marker),]
  print(table(temp$alleleA==whole$alleleA))
 
  assign(paste("prob", sp, sep="."), temp)
  print(table(temp$alleleA==whole$alleleA))
  
  ##Get difference
  diff <- (whole[,-(1:3)]-temp[,-(1:3)])

  map[,paste("sumdiff",sp,sep=".")] <- apply(abs(diff),1,sum)
  map[,paste("meandiff",sp,sep=".")] <- apply(abs(diff),1,mean)
  map[,paste("meddiff",sp,sep=".")] <- apply(abs(diff),1,median)
  map[,paste("per1diff",sp,sep=".")] <- apply(abs(diff)>.1,1,mean)
  map[,paste("per2diff",sp,sep=".")] <- apply(abs(diff)>.2,1,mean)
  map[,paste("per5diff",sp,sep=".")] <- apply(abs(diff)>.5,1,mean)
  
}

##Compare to random seed - baseline of random variation
table(whole$alleleA==whole2$alleleA)
table(names(whole)==names(whole2))
table(whole$marker==whole2$marker)
diff <- whole[,-(1:3)]-whole2[,-(1:3)]
sp <- "seed"
map[,paste("sumdiff",sp,sep=".")] <- apply(abs(diff),1,sum)
map[,paste("meandiff",sp,sep=".")] <- apply(abs(diff),1,mean)
map[,paste("meddiff",sp,sep=".")] <- apply(abs(diff),1,median)
map[,paste("per1diff",sp,sep=".")] <- apply(abs(diff)>.1,1,mean)
map[,paste("per2diff",sp,sep=".")] <- apply(abs(diff)>.2,1,mean)
map[,paste("per5diff",sp,sep=".")] <- apply(abs(diff)>.5,1,mean)




##Get breaks
segs <- gsub(".*[[:punct:]]", "", gsub("mb.*", "mb", gprob))
setwd(indir)
segdat <- data.frame("break"=1:3,stringsAsFactors=F)
for (seg in segs) {

  m1 <- read.table(paste(seg, "/1/1/markers_segment1.txt",sep=""), header=F, as.is=T)
  m2 <- read.table(paste(seg, "/1/2/markers_segment2.txt",sep=""), header=F, as.is=T)
  m3 <- read.table(paste(seg, "/1/3/markers_segment3.txt",sep=""), header=F, as.is=T)
  m4 <- read.table(paste(seg, "/1/4/markers_segment4.txt",sep=""), header=F, as.is=T)
  m5 <- read.table(paste(seg, "/1/5/markers_segment5.txt",sep=""), header=F, as.is=T)
  m6 <- read.table(paste(seg, "/1/6/markers_segment6.txt",sep=""), header=F, as.is=T)
  m7 <- read.table(paste(seg, "/1/7/markers_segment7.txt",sep=""), header=F, as.is=T)
  m8 <- read.table(paste(seg, "/1/8/markers_segment8.txt",sep=""), header=F, as.is=T)
  m9 <- read.table(paste(seg, "/1/9/markers_segment9.txt",sep=""), header=F, as.is=T)
  m10 <- read.table(paste(seg, "/1/10/markers_segment10.txt",sep=""), header=F, as.is=T)
  m11 <- read.table(paste(seg, "/1/11/markers_segment11.txt",sep=""), header=F, as.is=T)
  m12 <- read.table(paste(seg, "/1/12/markers_segment12.txt",sep=""), header=F, as.is=T)
  m13 <- read.table(paste(seg, "/1/13/markers_segment13.txt",sep=""), header=F, as.is=T)
  m14 <- read.table(paste(seg, "/1/14/markers_segment14.txt",sep=""), header=F, as.is=T)
	
  m1$seg <- 1
  m2$seg <- 2
  m3$seg <- 3
  m4$seg <- 4
  m5$seg <- 5
  m6$seg <- 6
  m7$seg <- 7
  m8$seg <- 8
  m9$seg <- 9
  m10$seg <- 10
  m11$seg <- 11
  m12$seg <- 12
  m13$seg <- 13
  m14$seg <- 14
	
  m <- rbind(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14)
  m$bp <- map$bp[match(m$V1, map$rsid)]
  m <- m[!is.na(m$bp),]
  minseg <- tapply(m$bp, m$seg, min)
  maxseg <- tapply(m$bp, m$seg, max)
  midseg <- (maxseg[-14]+minseg[-1])/2
  segdat[,paste("seg", seg, sep=".")] <- midseg
}

outdir <- "/data1/bsi/BORA_processing/devel/data_HTHGU/illumina1/MayoGAP/mayogap_emerge/raw/clean/beagle/run_beagle/chr22/result_new_paper_graph_11_01_2011"
setwd(outdir)
map <- map[!is.na(map$bp),]
png(file="chr10_median_dosage_diff_allsamples.png",res=100, width=1000, height=800)
#yrange <- range(c(map$meddiff.1mb, map$meddiff.2mb, map$meddiff.halfmb, map$meddiff.seed,),na.rm=T)
yrange <- c(0,0.20)
plot(lowess(map$bp, map$meddiff.seed,f=0.00001), type='l', ylim=yrange, main="Median of Absolute Value of Dosage Difference", xlab="BP", ylab="Median Absolute Difference")
lines(lowess(map$bp, map$meddiff.halfmb,f=0.00001),lty=4, col="green")
lines(lowess(map$bp, map$meddiff.1mb,f=0.00001),lty=2, col="red")
lines(lowess(map$bp, map$meddiff.2mb,f=0.00001), lty=3, col="blue")
lines(lowess(map$bp, map$meddiff.quatermb,f=0.00001), lty=5, col="cyan")
lines(lowess(map$bp, map$meddiff.halfquatermb,f=0.00001), lty=6, col="purple")
lines(lowess(map$bp, map$meddiff.quaterquatermb,f=0.00001), lty=7, col="magenta")

lines(lowess(map$bp[!is.na(map$r2)], yrange[2]-yrange[2]*map$r2[!is.na(map$r2)],f=0.00001), lty=5, col="orange")
abline(v=segdat$seg.1mb, col="red", lty=2)
abline(v=segdat$seg.2mb, col="blue", lty=3)
abline(v=segdat$seg.halfmb, col="green", lty=4)
abline(v=segdat$seg.quatermb, col="cyan", lty=5)
abline(v=segdat$seg.halfquatermb, col="purple", lty=6)
abline(v=segdat$seg.quaterquatermb, col="magenta", lty=7)

legend("topleft", col=c("black", "green", "red", "blue","cyan","purple","magenta"), lty=c(1,4,2,3,5,6,7), c("Different Seed", "Half MB", "1 MB", "2 MB","QUAT MB","HF_QUAT MB","QU_QUAT MB"), bty='n')
dev.off()
